#*********************************************************#
#                                                         #
#                                                         #
#                 REPLICATION FILES FOR:                  #
#                                                         #
#                 Barcelo and Labzina                     #
#                                                         #
#   Islamic State's Terrorist Attacks Disengage their     #
#     Supporters: Robust Evidence from Twitter            #
#                                                         #
#         British Journal of Political Science            #
#                                                         #
#                                                         #         
#*********************************************************#


#***************#
# Load Packages #
#***************#

library(foreign) 
library(boot)
library(stargazer)
library(MASS)
library(DirectEffects)
library(dplyr)

#***************#
#   Load Data   #
#***************#

#Set Working Directory

path <- ""

setwd(path)

#Load original BL datasets

re_data <- read.csv2("re_data.csv", sep = ",")
fe_data <- read.csv2("fe_data.csv", sep = ",")

re_data[,6:16] <- sapply(sapply(re_data[,6:16],as.character), as.numeric)
fe_data[,6:16] <- sapply(sapply(fe_data[,6:16],as.character), as.numeric)

################################
#Table 1: Omitted Variable Bias#
################################

main.sensitivity <- sensemakr(model = ate_mod50, 
                              treatment = "dm_Deaths_EurUS50_00",
                              benchmark_covariates = c("dm_nsuspensions_000"),
                              kd = c(1:11),
                              ky = c(1:11))

main.sensitivity_nonEur <- sensemakr(model = ate_mod50, 
                                     treatment = "dm_Deaths_nonEurUS50_00",
                                     benchmark_covariates = c("dm_nsuspensions_000"),
                                     kd = c(1:11),
                                     ky = c(1:11))

summary(main.sensitivity)
summary(main.sensitivity_nonEur)

###############################################################
#Table 2: ACDE of Terrorist Attacks on the Number of Followers#
###############################################################

fe_data$time_cat <- ifelse(fe_data$time > 122, 123, fe_data$time)
fe_data$time_cat <- as.factor(fe_data$time_cat)
fe_data$time_cat <- relevel(fe_data$time_cat, '1')

fe_data_deamean <- fe_data %>%
  group_by(account_id) %>%
  mutate(MeanValue50 = mean(Deaths_EurUS50_00), 
         dm_Deaths_EurUS50_00 = Deaths_EurUS50_00 - MeanValue50,
         MeanValue50_nonE = mean(Deaths_nonEurUS50_00), 
         dm_Deaths_nonEurUS50_00 = Deaths_nonEurUS50_00 - MeanValue50_nonE,
         MeanValue75 = mean(Deaths_EurUS75_00), 
         dm_Deaths_EurUS75_00 = Deaths_EurUS75_00 - MeanValue75,
         MeanValue75_nonE = mean(Deaths_nonEurUS75_00), 
         dm_Deaths_nonEurUS75_00 = Deaths_nonEurUS75_00 - MeanValue75_nonE,
         MeanValue100 = mean(Deaths_EurUS100_00), 
         dm_Deaths_EurUS100_00 = Deaths_EurUS100_00 - MeanValue100,
         MeanValue100_nonE = mean(Deaths_nonEurUS100_00), 
         dm_Deaths_nonEurUS100_00 = Deaths_nonEurUS100_00 - MeanValue100_nonE,
         MeanValue_foll = mean(lnfollowers_count),
         dm_lnfollowers_count = lnfollowers_count - MeanValue_foll,
         MeanValue_susp = mean(nsuspensions_000),
         dm_nsuspensions_000 = nsuspensions_000 - MeanValue_susp,
         MeanValue_rep = mean(Nreports_0000),
         dm_Nreports_0000 = Nreports_0000 - MeanValue_rep
         )

#Direct effects of the models using 50% discount rate

summary(ate_mod50 <- lm(dm_lnfollowers_count ~ dm_Deaths_EurUS50_00 + dm_Deaths_nonEurUS50_00 + 
             dm_nsuspensions_000 + dm_Nreports_0000 + 
             as.factor(time_cat) , data = fe_data_deamean))

form_main50 <- dm_lnfollowers_count ~ dm_Deaths_EurUS50_00 + dm_Deaths_nonEurUS50_00 + 
  as.factor(time_cat) | 1 |
  dm_nsuspensions_000 + dm_Nreports_0000

direct50 <- sequential_g(form_main50, data = fe_data_deamean)

summary(direct50)

#Direct effects of the models using 75% discount rate

summary(ate_mod75 <- lm(dm_lnfollowers_count ~ dm_Deaths_EurUS75_00 + dm_Deaths_nonEurUS75_00 + 
                          dm_nsuspensions_000 + dm_Nreports_0000 + 
                          as.factor(time_cat) , data = fe_data_deamean))

form_main75 <- dm_lnfollowers_count ~ dm_Deaths_EurUS75_00 + dm_Deaths_nonEurUS75_00 + 
  as.factor(time_cat) | 1 |
  dm_nsuspensions_000 + dm_Nreports_0000

direct75 <- sequential_g(form_main75, data = fe_data_deamean)

summary(direct75)

#Direct effects of the models using 100% discount rate

summary(ate_mod100 <- lm(dm_lnfollowers_count ~ dm_Deaths_EurUS100_00 + dm_Deaths_nonEurUS100_00 + 
                           dm_nsuspensions_000 + dm_Nreports_0000 + 
                           as.factor(time_cat) , data = fe_data_deamean))

form_main100 <- dm_lnfollowers_count ~ dm_Deaths_EurUS100_00 + dm_Deaths_nonEurUS100_00 + 
  as.factor(time_cat) | 1 |
  dm_nsuspensions_000 + dm_Nreports_0000

direct100 <- sequential_g(form_main100, data = fe_data_deamean)

summary(direct100)

summary(ate_mod100_sens <- lm(dm_lnfollowers_count ~ dm_Deaths_EurUS100_00 + dm_Deaths_nonEurUS100_00 + 
                           dm_nsuspensions_000 + 
                           as.factor(time_cat) , data = fe_data_deamean))

form_main100_sens <- dm_lnfollowers_count ~ dm_Deaths_EurUS100_00 + dm_Deaths_nonEurUS100_00 + 
  as.factor(time_cat) | 1 |
  dm_nsuspensions_000

direct100_sens <- sequential_g(form_main100_sens, data = fe_data_deamean)

summary(direct100_sens)

#We expand the g-sequential model to run sensitivity analysis on its main assumption using the discount rate of 100%

out_sens100 <- cdesens(direct100_sens, var = "dm_Deaths_EurUS100_00",
                       rho = c(cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*1,#1 times the observed correlation is -0.007
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*2,#2 times the observed correlation is -0.013
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*3,#3 times the observed correlation is -0.019
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*4,#4 times the observed correlation is -0.026
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*5,#5 times the observed correlation is -0.033
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*6,#6 times the observed correlation is -0.039
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*7,#7 times the observed correlation is -0.046
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*8,#8 times the observed correlation is -0.052
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*9,#9 times the observed correlation is -0.058
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*10,#10 times the observed correlation is -0.065
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*11,#11 times the observed correlation is -0.072
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*12,#12 times the observed correlation is -0.079
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*13,#13 times the observed correlation is -0.085
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*14,#14 times the observed correlation is -0.092
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*15,#15 times the observed correlation is -0.098
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*16,#16 times the observed correlation is -0.105
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*17,#17 times the observed correlation is -0.111
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*18,#18 times the observed correlation is -0.118
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*19,#19 times the observed correlation is -0.124
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*20,#20 times the observed correlation is -0.131
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*21,#21 times the observed correlation is -0.138
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*22,#22 times the observed correlation is -0.144
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*23,#23 times the observed correlation is -0.151
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*24,#24 times the observed correlation is -0.157
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*25,#25 times the observed correlation is -0.164
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*26,#26 times the observed correlation is -0.170
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*27,#27 times the observed correlation is -0.177
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*28,#28 times the observed correlation is -0.183
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*29,#29 times the observed correlation is -0.19
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*30#30 times the observed correlation is -0.197
                                                      ), 
                       bootstrap = "standard", 
                       boots_n = 100, verbose = TRUE)

plot(out_sens100)

plot_outsens100 <- out_sens100

plot_outsens100_data <- as.data.frame(cbind(out_sens100$rho, out_sens100$acde, out_sens100$se))
colnames(plot_outsens100_data) <- c('rho', "acde", "se")

########################################################################
# Figure 1: Sensitivity Analysis of the ACDE of Attacks on Followers
#     Figure 1A: Europe and the US attacks
#     Figure 1B: Attacks elsewhere
########################################################################

acde_sens <- ggplot(plot_outsens100_data,
       aes(x = rho,
           y = acde)) +
  # Add a ribbon with the confidence band
  geom_smooth(colour = "black",
              aes(
                # lower and upper bound of the ribbon
                ymin = acde-1.96*se, 
                ymax = acde+1.96*se
              ),
              stat = "identity") +
  geom_hline(
    yintercept = 0,
    linetype='dotted'
  ) +
  xlab("Correlation between mediator and outcome errors") +
  ylab("Estimated ACDE of the Number of Victims in Europe and the U.S.") + 
  scale_x_continuous(
    breaks = c(plot_outsens100_data$rho[1], plot_outsens100_data$rho[5], 
               plot_outsens100_data$rho[10], plot_outsens100_data$rho[15],
               plot_outsens100_data$rho[20], plot_outsens100_data$rho[25],
               plot_outsens100_data$rho[30]),
    
    labels=c('-0.20' = '-0.20 \n = 30 X cor(M1,Y)',
             '-0.16' = '-0.16 \n = 25 X cor(M1,Y)', 
             '-0.13' = '-0.13 \n = 20 X cor(M1,Y)', 
             '-0.10' = '-0.10 \n = 15 X cor(M1,Y)', 
             '-0.07' = '-0.07 \n = 10 X cor(M1,Y)', 
             '-0.03' = '-0.03 \n = 5 X cor(M1,Y)', 
             '-0.01' = '0.00 \n = 1 X cor(M1,Y)'
    )) +
  theme_minimal()

png(filename="acde_sens.png", 
    units="in", 
    width=9, 
    height=6, 
    pointsize=15, 
    res=1000)
acde_sens
dev.off()

out_sens100n <- cdesens(direct100_sens, var = "dm_Deaths_nonEurUS100_00",
                       rho = c(cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*1,#1 times the observed correlation is -0.007
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*2,#2 times the observed correlation is -0.013
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*3,#3 times the observed correlation is -0.019
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*4,#4 times the observed correlation is -0.026
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*5,#5 times the observed correlation is -0.033
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*6,#6 times the observed correlation is -0.039
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*7,#7 times the observed correlation is -0.046
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*8,#8 times the observed correlation is -0.052
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*9,#9 times the observed correlation is -0.058
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*10,#10 times the observed correlation is -0.065
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*11,#11 times the observed correlation is -0.072
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*12,#12 times the observed correlation is -0.079
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*13,#13 times the observed correlation is -0.085
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*14,#14 times the observed correlation is -0.092
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*15,#15 times the observed correlation is -0.098
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*16,#16 times the observed correlation is -0.105
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*17,#17 times the observed correlation is -0.111
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*18,#18 times the observed correlation is -0.118
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*19,#19 times the observed correlation is -0.124
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*20,#20 times the observed correlation is -0.131
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*21,#21 times the observed correlation is -0.138
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*22,#22 times the observed correlation is -0.144
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*23,#23 times the observed correlation is -0.151
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*24,#24 times the observed correlation is -0.157
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*25,#25 times the observed correlation is -0.164
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*26,#26 times the observed correlation is -0.170
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*27,#27 times the observed correlation is -0.177
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*28,#28 times the observed correlation is -0.183
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*29,#29 times the observed correlation is -0.19
                               cor(fe_data_deamean$lnfollowers_count, fe_data_deamean$nsuspensions_000)*30#30 times the observed correlation is -0.197
                       ), 
                       bootstrap = "standard", 
                       boots_n = 100, verbose = TRUE)

plot(out_sens100n)

plot_outsens100n <- out_sens100n

plot_outsens100n_data <- as.data.frame(cbind(out_sens100n$rho, out_sens100n$acde, out_sens100n$se))
colnames(plot_outsens100n_data) <- c('rho', "acde", "se")


acde_sens_nonEur <- ggplot(plot_outsens100n_data,
       aes(x = rho,
           y = acde)) +
  # Add a ribbon with the confidence band
  geom_smooth(colour = "black",
              aes(
                # lower and upper bound of the ribbon
                ymin = acde-1.96*se, 
                ymax = acde+1.96*se
              ),
              stat = "identity") +
  geom_hline(
    yintercept = 0,
    linetype='dotted'
  ) +
  xlab("Correlation between mediator and outcome errors") +
  ylab("Estimated ACDE of the Number of Victims outside Europe and the U.S.") + 
  scale_x_continuous(
    breaks = c(plot_outsens100n_data$rho[1], plot_outsens100n_data$rho[5], 
               plot_outsens100n_data$rho[10], plot_outsens100n_data$rho[15],
               plot_outsens100n_data$rho[20], plot_outsens100n_data$rho[25],
               plot_outsens100n_data$rho[30]),
    
    labels=c('-0.20' = '-0.20 \n = 30 X cor(M1,Y)',
             '-0.16' = '-0.16 \n = 25 X cor(M1,Y)', 
             '-0.13' = '-0.13 \n = 20 X cor(M1,Y)', 
             '-0.10' = '-0.10 \n = 15 X cor(M1,Y)', 
             '-0.07' = '-0.07 \n = 10 X cor(M1,Y)', 
             '-0.03' = '-0.03 \n = 5 X cor(M1,Y)', 
             '-0.01' = '0.00 \n = 1 X cor(M1,Y)'
    )) +
  theme_minimal()

png(filename="acde_sens_nonEur.png", 
    units="in", 
    width=9, 
    height=6, 
    pointsize=15, 
    res=1000)
acde_sens_nonEur
dev.off()

png(filename="acde_sens_all.png", 
    units="in", 
    width=10, 
    height=12, 
    pointsize=15, 
    res=1000)
ggarrange(acde_sens, acde_sens_nonEur,
          labels = c("A) ", "B) "),
          ncol = 1, nrow = 2)
dev.off()

